De novo mitochondrial genome sequencing of Cladonia subulata and phylogenetic analysis with other dissimilar species

In this study, the complete mitochondrial genome of Cladonia subulata (L.) FH Wigg was sequenced and assembled and then compared with those of other Cladonia species. The mitogenome of Cladonia subulata, the type species of Cladonia, consisted of a circular DNA molecule of 58,895 bp 44 genes (15 protein-coding genes, 2 rRNA genes, and 27 tRNA genes). The base composition had shown an obvious AT preference, and all 27 tRNA genes formed a typical clover structure. Comparison with other 7 Cladonia species indicated that the duplication/loss of tRNAs had occurred during evolution, and introns appeared to explain the variation in cox1 genes in Cladonia, the mitochondrial genome tends to be generally conservative and local dynamic changes. Repeat sequences were mainly located in gene intervals, which were mainly distributed among intergenic spacers and may cause rearrangement of the mitogenome. The phylogenetic results showed that Cladonia subulata and C. polycarpoides were assigned to the Cladonia Subclade. The results add to the available mitochondrial genome sequence information of Cladonia subulata, provide basic data for the systematic development, resource protection, and genetic diversity research in Cladonia subulata, and also provide theoretical support for further genomic research of lichens.


Introduction
Mitochondria play an important role in eukaryotes by producing adenosine triphosphate (ATP) through oxidative phosphorylation, and mitochondria contain their own genome. In contrast to the nuclear genome, the mitochondrial genome usually had its own coding system and follows the rules of uniparental inheritance, which is not affected by the recombination of two parental chromosomes [1]. Because of its rapid evolution and high copy number per cell, mitogenome DNA has been widely used to study the origin, evolution, and phylogeny of fungi [2][3][4]. The fungal mitogenome DNA shows a high DNA AT-content and a wide range of genome sizes [5]. The fungal mitogenome contains 14 protein-coding genes (PCGs) (nad1-6, nad4L, cob, cox1-3, atp6, atp8 and atp9) that are related to electron transport and oxidative phosphorylation in the inner membrane of mitochondria. It also includes the mitochondrial 16S ribosomal small subunit RNA gene (rrns), the mitochondrial 23S ribosomal large subunit RNA gene (rrnl) and a variety of tRNA genes necessary for translation [6]. Lichens, are pioneer colonizers, that play a major role in studies aimed reconstructing early life on land [7][8][9][10]. A recent study suggests that lichens are self-sustaining ecosystems main consisting of Ascomycota or Basidiomycota fungi and prokaryotic or eukaryotic algae, and other microscopic organisms [11,12]. This special symbiotic system allows lichens to survive in extreme environments such as waterless deserts or areas of extreme cold [13]. In addition, lichens also play an extremely important role in antibacterial and antitumor activities because of their rich lichenic acid [14][15][16][17][18][19].
Lichens show rich morphological and ecological diversity in nature [20], however, the evidence based on morphological and chemical characteristics is not enough to accurately reflect the species diversity or species delimitation, and may misrepresent the diversity of lichenized fungi, especially in closely related species that have very similar morphological, the combination of molecular data can more accurately classify species [21][22][23]. Mitogenome is not only conducive to solving the problems of taxonomic identification, especially some symbiotes (such as arbuscular mycorrhizal fungi) that are difficult to be characterized by conventional methods [24], but also conducive to the comparative analysis of mitogenome to reveal the evolutionary events between species, such as the dynamic changes of genes [25].
Cladonia is one of the largest genera of lichens with approximately 475 species showing a variety of morphologies and habits [26], and this genus has been used in ecological research as well as in research related to industry and medicine [27][28][29][30]. To date (2023.02.01), almost 17 complete mitogenomes of this genus have been released in the NCBI database, among which the complete mitogenome sizes from 46.2 to 66.1 kb, and few of their mitogenes have been annotated or analysed, the available mitogenome data for Cladonia are far less abundant than those for nuclear genomes (https://www.ncbi.nlm.nih.gov/nuccore), so there are still some gaps in Cladonia phylogenetics and comparative genomics.
Cladonia subulata is the type species of Cladonia [31], and has been included in the "Chinese Biodiversity Red List-Macrofungi" published by the Chinese Academy of Sciences, Ministry of Ecology and Environment, as a species of LC grade (www.mee.gov.cn). In a recent study, it was found that Cladonia subulata showed strong cytotoxic activity which could provide a potential treatment for cancer [19]. Cladonia subulata has many morphological characters and can be distinguished based on lichen acids, podetia, and internal transcribed spacer sequences (ITS) [32]. Compared with the simple nuclear genome fragment, the complete mitochondrial genome could provide more abundant genetic resources for these protected species and conducive to the exploration of its functional genes [33]. The composition, repeat sequence, gene arrangement and other differences in mitogenomes could reflect the genetic evolution level in species [34], which is beneficial to the protection of the diversity of biological genetic resources and the study of species population structure by mitogenomes [24,35]. However, all of nucleotide records of Cladonia subulata in NCBI commonly used in phylogenetic and/or molecular systematic studies, the research on the mitogenome of C. subulata is still lacking. This information gap limits our understanding of Cladonia subulata at a genetic level.
In this study, we first assembled and annotated the complete mitogenome of Cladonia subulata. To reveal the characteristics of its mitogenome, and its differences from those of other species in the genus, the mitogenomes of Cladonia subulata and seven other species in the genus were compared.
Republic of China on the protection of wildlife'. The voucher specimen Cladonia subulata (YH0097) was deposited in the "Lichen's Research Center in Arid Zones of Northwest China", Xinjiang University, Xinjiang, China. We first identified this specimen based on morphology [20,32]: Squamules greyish green to green, esorediate or sometimes whitish sorediate at the margin. Potedia abundant, 15-40 mm tall, usually unbranched, greyish green, subulate but finally forming cup at the tip. Surface ecorticate, covered by powdery soredia, and small podetial squamules occurring at the base. Apothecia seldom, stalked, brown. Spot tests: K-, C-, KC-, P +red. Secondary metabolite: fumarprotocetraric acid (S1 Fig).
Total genomic DNA was extracted using the Fungi Genomic DNA Extraction Kit (Solarbio, Beijing, China) and stored at -20˚C. Then, rDNA-ITS, mtSSU and RPB2 were selected as the gene markers via PCR, rDNA-ITS was used a pair of primers ITS1F [36] and ITS4 [37]; for mtSSU the primer pair mtSSUl and mtSSU3R [38]; the RPB2 was amplified using nested PCR [34], the primers RPB2-5F and RPB2-7R [39] were used firstly, RPB2dRaq and RPB2rRaq [32] was used for second reaction. And their conditions of PCR were following the S1 Table. Then, the PCR products were sequenced by Sangon Biotech Co. Ltd. (Shanghai, China). The ITS, mtSSU and RPB2 sequences were submitted to NCBI under accession numbers ON920702, ON937425 and OQ025089. After sample QC, the gDNA was used to construct a singlestranded circular (ssCir) library, and the ssCir library was then amplified through rolling circle amplification (RCA) to obtain DNA nanoballs (DNBs), which were loaded into flow cells and sequenced on the DNBSEQ Platform by BGI Biotechnology Co. (Shenzhen, China). Approximately 6 Gb of raw data were generated for each library, and the total read Q30 was over 80%.

Phylogenetic analysis
To further prove from a molecular point of view that the species was Cladonia subulata. The ITS, mtSSU and RPB2 sequences were concatenated to verify C. subulata, Pseudevernia cladonia as an outgroup and three loci in Cladonia were extracted from GeneBank (S2 Table) followed by BLASTN with corresponding C. subulata genes using MAFFT [54], Gblocks was used to remove ambiguous regions and keep conserved regions of the sequences. Best-fit model according to Bayesian Information Criterion (BIC) by ModelFinder [55] which were used to the MrBayes: TIM3e+G4: ITS+RPB2; K3Pu+F+G4: SSU. Bayesian analyses were performed using MrBayes v.3.1.2 [56] with the following main parameters: ngen = 1000000, nchains = 4, simple freq = 100, nst = 6, rates = gamma, burn in = 0.25. The tree was viewed in The complete mitogenomes of Cladonia subulata and 23 Lecanorales species from Gen-Bank (S3 Table) were used for phylogenetic comparison. Heterodermia casarettiana (A. Massal.) Trevis. (NC 042185) and H. speciosa (Wulfen) Trevis. (NC 040159) were used as outgroup taxa [58,59]. The phylogenetic analysis used 15 PCGs of mitogenomes, which were first aligned using MAFFT. After editing and trimming to produce a sequence matrix, the results were concatenated by PhyloSuite [60]. ModelFinder [56] was used to select the best-fit partition model of each PCG (edge-linked) for maximum likelihood (ML) and Bayesian analyses by using the AIC (S4 Table). An ML analysis was performed using IQ-tree [61] incorporated in PhyloSuite under the Edge-linked partition model with 5000 ultrafast bootstraps replicates. Subsequently, Bayesian analyses were conducted using the Markov chain Monte Carlo (MCMC) method in MrBayes v.3.1.2 according to the Kundu [62] with the following conditions: the partition model was chosen (S4 Table), then the tree was then run for 10,000,000 generations with 25% burn in, with trees saved saving at every 100 generations. MCMC analysis was used to generate the convergence metrics until the standard deviation (SD) of split frequencies fell below 0.01 and the potential scale reduction factor (PSRF) for all parameters approached 1.0. Both BI and ML phylogenetic trees were viewed and edited via the web based iTOL tool (https://itol.embl.de/) [63].

Mitogenomic characterization of Cladonia subulata
The complete mitogenome of Cladonia subulata was assembled as a circular molecule of 58,895 bp (Fig 1), including 15 PCGs, seven for NAD dehydrogenases (nad1, nad2, nad3, nad4, nad4L, nad5, nad6), three for cytochrome oxidases (cox1, cox2, cox3), three for ATP synthases (atp6, atp8, atp9), one for cytochrome b (cob) and one for ribosomal protein subunit 3 (rps3), 2 ribosomal RNA genes (rrnS, rrnL) and 27 transfer RNA genes ( Table 1). All genes were transcribed at the same strand. The length of 15 PCGs varied from 147 bp (atp8) to 8,208 (cox1) and the total length was 26,194 bp which accounted for 44.48% of the mitogenome. The length of tRNAs was 1989 bp, being 3.38% of the mitogenome, while rRNAs was 2,420 bp, which was 4.11% of the mitogenome. Previous study [64] has indicated that the intron could influence the size of lichenized fungi. In our study, the mitogenome contained 11,583 bp introns in PCGs, and the introns accounted for 44.22% of the PCGs and 19.67% of the mitogenome, especially in cox1, which contained 79.57% introns.
The GC content of mitogenomes varies among different species, which may be affected by mutation bias, selection and recombination-related DNA repair bias [65]. The nucleotide composition of Cladonia subulata showed an AT bias (70.70%) and the GC content was 29.30%, similar to most Cladonia species (Tables 1 and 2). According to the second parity rule of Chargaff (A-T and G-C within one stand) [66], the mitogenome of Cladonia subulata presented a negative AT skew (-0.015) and a positive GC skew (0.068) ( Table 1). This indicates that the mitogenome of Cladonia subulata contains more T than A nucleotides and more G than C nucleotides, as reported in other mitogenomes of Cladonia.
Four PCGs (nad4, rps3, cob, and cox2) were initiated with GTG, TTA, ATA, and GTG respectively. Other PCGs were started with the typical ATG start codons. The majority of PCGs were terminated by typical TAA codons, except for rps3, cob, atp6 and nad1 which had terminal codons TAG (Table 1). In addition, the first nucleotide (A) of the nad5 initiation codon ATG was also the last nucleotide of the nad4L termination codon TAA (Tables 1 and  S5). All protein and RNA genes were encoded on the positive strand, the trnL was included in cox1 (Fig 1, S5 Table). All rRNAs were independent, and the longest distance was 1,581 (rrnS). Subsequently, 16 tRNAs located both sides of rrnS (S5 Table). This phenomenon suggests that

PLOS ONE
the mitogenome of Cladonia subulata has a non-continuous segment and lack of overlapping genes (S5 Table). tRNAs present the same significance during translation as mRNAs and proteins [67]. In Cladonia subulata, 27 tRNA genes in the mitogenome encoded the 20 standard amino acids, and the secondary structures of all tRNAs were successfully predicted ranging from 69 bp to 86 bp (Fig 2, S5 Table). Different types of analyses may result in differences in the predicted

PLOS ONE
typical tRNA secondary structure. In addition, 27 tRNAs exhibited typical clover-leaf secondary structure; e.g., trnY (Tyr), trnS (Ser), and trnL (Leu) showed a variable loop. All tRNAs presented typical clover-leaf secondary structures, and the remaining mismatched bases had showed non-canonical G-U pairs (Fig 2). In accordance with the analysis of secondary structure in fungi, the variations in the numbers of extra arms may cause differences in tRNA length.
Codons are carriers used for identifying and transmitting the genetic information of organisms and play an important role in biological genetics and variation [68]. The codon usage patterns of fungal mitochondria are generally similar, including the most commonly used codons UUU, UUA, AUU, AUG, GUU, CUU, AUA, and GUA, and 30 other codons, with exception of UAA. According to the analysis of codon usage, the most frequently used codons were AGA (for Arg), UUA (for Leu), CCU (for Pro). It is observed that the frequency of A and U in the codon is very high (Fig 3, S6 Table).

Comparative mitogenomes of Cladonia fungi
The seven published and annotated fungal mitogenomes of Cladonia were compared to determine the variation in fungal mitogenomes. In these 8 species, the mitogenome size from 45,312 to 59,878 bp. The length of PCGs was the main reason for the differences in mitogenome length ( Table 2). The tRNAs number various from 19 to 27, which indicated that some species may have lost genes in the course of evolution. The loss of mt-tRNA genes may be associated with selection pressure or transfer to the nucleus [69]. The numbers and sizes of introns in the 8 species of Cladonia were measured and compared, the size of introns from 2,617 to 13,642 which showed a big difference. In these species, the largest values were obtained in Cladonia rangiferina which contained 17 introns with a total length of 13,642 bp while the smallest values were obtained in Cladonia peziziformis which contained 4 introns with a total size of 2,617 bp.

PLOS ONE
The synteny results obtained from in Mauve showed that the homologous collinear blocks occurred in all 8 fungal mitogenomes of Cladonia (Fig 4). The PCGs and rRNAs of the 8 species were obviously highly similar in composition. The 15 PCGs of these 8 species were cox1 -nad1 -nad4 -rps3 -cob -cox2 -atp9 -atp6 -atp8 -nad6 -cox3 -nad2 -nad3 -nad4L -nad5, and all genes were arranged in the same order in these 8 Cladonia species (S7 Table). At the synteny level, the mitogenome of Cladonia subulata showed a conserved in the gene arrangement and gene order of 15 PCGs, 2 rRNA genes and the majority of tRNA genes relative to those generally observed in Cladonia. Subsequently, in the mitochondrial genome of Cladonia subulata, it is divided into 6 gene clusters (Fig 4). The locations of the gene clusters were basically the same, but there were some differences in size. This main difference was located near cox1, which may be due to the length differences among homologous clusters caused by the acquisition or deletion of introns, resulting in a change in mitochondrial genome length (Fig 4). To better understand whether the change in homologous clusters was highly variable in litmus, 18 public complete mitochondrial genome sequences (including unannotated sequences) of Cladonia in NCBI were selected for further collinearity analysis. A total of 13 homologues were clearly identified the 18 mitochondrial genomes (S3 Fig). It was found not only that the mitochondrial genome of Cladonia showed a rearrangement phenomenon, but also that some homologous clusters were lost when more diverse samples were analysed. Interestingly, only Cladonia uncialis showed notable rearrangement and loss of homologous clusters, with four inverted and two lost homologous clusters. In a previous study, Cladonia uncialis showed inversion blocks of four genes (nad6, cox3, nad2 and cox3) compared with C. rangiferina [26]; however, combined with the results of gene arrangement and collinearity, there was no largescale rearrangement in the other 17 Cladonia species which may indicate that the mitogenome of Cladonia is relatively conservated.
A more in-depth comparative study of these 8 species is shown in Fig 5. Among the 15 PCGs, the lengths of atp8, atp9, cox3, and nad3 remained the same, while the lengths of other genes showed varying degrees of divergence (Fig 5A). The gene with the greatest length difference was cox1, which was related to the abundant introns it contained. The GC content in the
Most PCGs had a positive G-C skew, while atp8, cox2 and nad2 showed a negative G-C skew. In addition, the G-C skew of rps3 in Cladonia rangiferina was zero which meant that the G and C were present in equal amounts, and only rps3 of C. subulata exhibited a positive G-C skew (Fig 5D).
The comparative analysis of tRNAs in Cladonia showed that C. subulata had two copies of trnL (Leu) and trnS (Ser), three copies of trnM (Met) and trnR (Arg) (S7 Table). Different species presented different tRNA copy numbers. Cladonia peziziform had three copies of trnM, and two copies of trnR, trnS, and trnL; Cladonia macilenta had no copies; Cladonia leporina had four copies of trnT, and three copies of trnM, and two copies of trnR, trnS, and trnL; Cladonia apodocarpa had three copies of trnM, and two copies of trnR, trnL, and trnK; Cladonia petrophila had two copies of trnM, trnR, trnS, and trnK; Cladonia rangiferina had three copies of trnR and two copies of trnL; and Cladonia subtenuis had three copies of trnN and trnR, and two copies of trnS and trnL (S4 Fig, S7 Table). The varying copy numbers in these species indicated that Cladonia macilenta may have lost at least one copy of trnM, trnR, trnS, trnL, trnT, and trnK or that duplications of these tRNAs have occurred in other species. Subsequently,

PLOS ONE
Cladonia subulata may have lost one copy of trnK and three copies of trnT, or these two tRNAs may have been duplicated in other species. In general, sequences different in tRNA genes were connected with variation.

Repeat sequences
BiBiserv-REPuter can locate and identify forwards, reverse, complementary and palindromic repeat sequences [52]. The accumulation of repetitive sequences in the fungal mitogenome over time may lead to dynamic changes in genome structure, and in turn affect the rearrangement of the mitogenome [70]. In the mitogenome of Cladonia subulata, REPuter identified 67 forward types (totalling 2,640 bp), 22 palindromic repeats (totalling 1,456 bp) and 1 reverse repeat of only 30 bp (S8 Table). Tandem repeats were searched by using Tandem Repeats Finder, which identified 30 tandem repeats in Cladonia subulta, and the longest tandem sequence of C. subulata was 42 bp. Most tandem repeats were located in the intergenic spacer, and some were located in the cox1 gene (Fig 6, S9 Table). Simple sequence repeats (SSRs) are a valuable type of molecular marker that may show high variation within a given species and has

PLOS ONE
been used in population genetics and polymorphism research [71]. We analysed the occurrence and types of SSRs in the mitochondrial genome of stamens, and a total of 93 SSRs were identified (Tables 3 and S10). Most of these SSRs were composed of single-nucleotide and double-nucleotide repeats, which were found 36 times and 50 times, respectively, and trinucleotide repeats (2), tetranucleotide repeats (3), and hexanucleotide repeats (2) were also present at low frequencies. All single-nucleotide repeats consisted of A/T repeats. Similarly, all dinucleotide repeats consisted of AT/AT repeats (Tables 3 and S10). The results showed that the mitochondrial genome SSRs were mainly composed of short PolA or PolyT repeats and rarely contained tandem G or C repeats. It can be seen from Fig 6 that most of the repetitions were located in the gene intervals.

Phylogenetic analysis
The two phylogenetic trees (maximum likelihood tree and Bayesian inference tree) lead to the same topology with different supporting value. Phylogenetic analysis based on 23 mitochondrial genome PCGs of Lecanorales revealed division into four clades with high support (Fig 7). The phylogenetic analysis showed that there was a closer genetic relationship between Parmeliaceae and Lecanoraceae. All the Cladonia species were well clustered in the monophyletic branch of the Cladoniaceae. The results were similar to the previous multigene Lecanoromycetes phylogenetic analysis [58,59]; however, these researches also indicated that a certain portion of the phylogeny was less stable between Parmeliaceae and Lecanoraceae due to a higher concentration of taxa with only two (RNA-coding) genes, which suggests that more sequence loci were beneficial for constructing a stable phylogenetic tree.
Finally, to obtain a better-supported phylogenetic tree of Cladonia, phylogenetic analysis was carried out on the complete mitogenomes of 18 Cladonia subulata (including annotated and nonannotated mitogenomes), Usnea mutabilis as the outgroup; and the conditions for performing the phylogenetic analysis were the same as those applied for Lecanorales (Fig 8).
Our phylogenetic results were similar to the most recent multiple-locus (nuclear and protein loci) phylogeny of the genus Cladonia and some previous studies [73][74][75][76]. According to the results (Fig 8), Cladonia subulata formed a well-supported clade with C. polycarpoides, both species have "non-classic" cup-shaped. Cladonia furcata formed a high supported clade with C. peziziformis, which belonged to the subclade Ascyphiferae according to they contained the fumarprotocetraric acid and sparsely or moderately branched. Both our results and Brigham's [75] showed that the Cladonia stipitata and C. petrophila had a more related relationship than C. apodocarpa in subclade Apodocarpae, most of them with non-calcareous substrates, ascending squamules, rare podetia, fumarprotocetraric acid present. Interestingly, Lendemer's [74]

Fig 7. Phylogeny of Lecanorales based on the 23 available mitochondrial genomes based on the Bayesian inference (BI) and maximum likelihood (ML).
Numbers on branches indicate posterior probability (BI) and bootstrap support (ML). Cladonia is framed with a dotted blue box. The species and GenBank accession numbers of the mitogenomes used in the phylogenetic analysis are provided in S3 Table. https://doi.org/10.1371/journal.pone.0285818.g007

PLOS ONE
study showed that Cladonia stipitata were morphologically similar to C. petrophila but C. stipitata formed a well-supported monophyletic group which distinct from C. apodocarpa and C. petrophila in phylogenetic results, the difference between us requires more mitogenomes of these three species for further phylogenetic analysis to get a more accurate phylogenetic result. In addition, Cladonia robbinsii and C. pyxidata belong to the subclade Foliaceae and subclade Graciles respectively. The subclade Cladonia, subclade Ascyphiferae, subclade Apodocarpae, subclade Foliaceae and subclade Graciles belong to Clade Cladonia. Cladonia leporina and C. ravenlii belonged to the Clade Erythrocarpae along with C. coccifera and C. macilenta [33,76] according to their red hymenial disc. Cladonia subtenuis and C. rangiferina formed a well-supported clade which belonged to Clade Crustaceae, the species in this clade were widely distributed in arctic and boreal to temperate areas. Cladonia uncialis, C. caroliniana and C. squamosa belong the Clade Unciales, Clade Borya, Clade Perviae according to the Stenroos's [73] results. In addition, all clades within the tree had sufficient bootstrap support (Fig 8), it was showed the Clade Erythrocarpae had related relationship with Clade Pervia, and the Clade Unciales, Clade Borya and Clade Perviae showed a more related relationship than other Clades. Over all, the PCGs showed a great application in Cladonia phylogenetic analysis. The sequencing of mitochondrial genes has been widely performed in fungal studies [77,78]. In lichenized fungi, although high variability of ITS sequences was found in the studies, SSU, LSU and cox genes have been widely used in taxonomic studies of lichenized fungi and have shown greater stability when used in phylogenetic construction in previous studies [79][80][81][82].

Conclusion
To our knowledge, this study was the first in which the mitogenome of Cladonia subulata was sequenced and assembled, and its overall characteristics were annotated and analysed, its mitogenome main characteristics were similar to those observed in other species of the genus. The cox1 gene showed the largest variation among all PCGs because of abundant introns. The presence/absence of introns indicated that introns were related to the evolutionary dynamics of mitotic introns [6]. Subsequently, the duplication of the tRNAs of the C. subulata mitogenome may have resulted in a number of repeat/loss events during evolution.
Comparative analysis of mitogenome alignment and structure can be used to reveal evolutionary relationships among species [83]. In the gene collinearity analysis indicated dynamic changes in the mitogenome of Cladonia. Cladonia uncialis showed a larger change in homologous clusters than the other species, giving rise to the question of whether this phenomenon is common or specific to Cladonia. Further statistical analysis will be required to answer this question. The repeat sequences of Cladonia subulata mitogenome were mainly distributed in gene intervals (Fig 6), which are more likely to show variation than coding sequences. The vast majority of sequences were composed of A/T nucleotides, which was consistent with the high A/T content (70.70%) of the genome. Therefore, whether a high A/T content is related to high variation is worthy of further studying.
This study provides new insights into the genetics of Cladonia and comparative analysis of eight fungal mitogenomes provides an understanding of the mitochondrial evolution in the genus, which will provide data support for further study and also protect the diversity of Cladonia genes to some extent.